%% Read distance data
clear all
close all
cd 'PATH_TO_PROJECT/Data/GS/gstime_all'
f = 'distance_matrix_germany_stid.csv'
idd = readmatrix(f);

d0 = readmatrix('d0.csv');
d1 = readmatrix('d1.csv');
d2 = readmatrix('d2.csv');
d3 = readmatrix('d3.csv');
d4 = readmatrix('d4.csv');
d5 = readmatrix('d5.csv');
d6 = readmatrix('d6.csv');
d7 = readmatrix('d7.csv');
d8 = readmatrix('d8.csv');
distance = [d0 d1 d2 d3 d4 d5 d6 d7 d8];
clear d0 d1 d2 d3 d4 d5 d6 d7 d8
% distance = readmatrix(f);

distance=distance-diag(diag(distance));
m=length(idd(:,1)); % Number of stations
m_original=m;
min_distance=zeros(m,1);
for t=1:m
min_distance(t)=min(setdiff(distance(t,:),min(distance(t,:))));
end

cd 'PATH_TO_PROJECT/Data/Market/all'
for monopoly_threshold = [20]
	disp(monopoly_threshold)
	%%%%% Define monopolies as stations at more than 20 minutes from their rivals. 
	index=(min_distance(:,1)<=monopoly_threshold);
	%'Index monopolies'
	Mon=idd(index==0,:);
	%%%% Select only stations that are not monopolies
	distance2=distance(index>0,index>0);
	idd2=idd(index>0,1);
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	m=length(idd2(:,1)); % Number of stations that are not monopolists
	Y=squareform(distance2);
	Z = linkage(distance2,'average');
	I = inconsistent(Z,3);

	for prune_threshold = [80]
		disp(prune_threshold)
		c = cluster(Z,'cutoff',prctile(I(:,4),prune_threshold), 'depth',3);
		%%% Data to export
		data_export=[idd2 c];
		label={'idd2','cluster'};
		f_out = strcat('germany_', num2str(monopoly_threshold), '_', num2str(prune_threshold));
		disp(f_out)
		xlswrite(f_out,data_export,'cluster','A2');
	end
end
